2021 Fasting Cancer Project

Set options and working directory

Load libraries

library(Seurat)
library(cowplot)
library(umap)
library(dplyr)
library(Matrix)
library(readxl)
library(openxlsx)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(sctransform)
library(knitr)
library(kableExtra)
#library(biomaRt)
library(DESeq2)
library(escape)
library(dittoSeq)
library(GSEABase)
library(scater)
library(ComplexHeatmap)

Import marker sets

gene.lists <- read_xlsx("Munoz_Yilmaz_CellCycle_signatures.xlsx")
gene.lists.names <- colnames(gene.lists)
GOI.lists <- c()
for (i in gene.lists.names) {
  tmpList <- gene.lists %>% dplyr::select(all_of(i))
  tmpList <- tmpList[!is.na(tmpList)]
  GOI.lists[[i]] <- tmpList
}

Load the Cell Ranger Matrix Data and create the base Seurat object.*

the initial processing was done with r 3.6.1 with Seurat 3.2.0 - the UMAP comes out slightly differently in r 4.0.3 with Seurat 3.2.3*

#al.dat <- Read10X("200218Yil_data/al/filtered_feature_bc_matrix/")
#f.dat <- Read10X("200218Yil_data/f/filtered_feature_bc_matrix/")
#rf.dat <- Read10X("200218Yil_data/rf/filtered_feature_bc_matrix/")
#rf.rapa.dat <- Read10X("200218Yil_data/rf.rapa/filtered_feature_bc_matrix/")

#al <- CreateSeuratObject(counts = al.dat, project = "al", min.cells = 3, min.features = 200)
#f <- CreateSeuratObject(counts = f.dat, project = "f", min.cells = 3, min.features = 200)
#rf <- CreateSeuratObject(counts = rf.dat, project = "rf", min.cells = 3, min.features = 200)
#rf.rapa <- CreateSeuratObject(counts = rf.rapa.dat, project = "rf.rapa", min.cells = 3, min.features = 200)

#al[["percent.mt"]] <- PercentageFeatureSet(al, pattern = "^mt-")
#f[["percent.mt"]] <- PercentageFeatureSet(f, pattern = "^mt-")
#rf[["percent.mt"]] <- PercentageFeatureSet(rf, pattern = "^mt-")
#rf.rapa[["percent.mt"]] <- PercentageFeatureSet(rf.rapa, pattern = "^mt-")

#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(rf, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

Merge samples to single Seurat Object

#merged_Seu <- merge(al, c(f,rf,rf.rapa), project = "Diet")

#merged_Seu <- NormalizeData(merged_Seu, normalization.method = "LogNormalize", scale.factor = 10000)
#merged_Seu <- FindVariableFeatures(merged_Seu, selection.method = "vst", nfeatures = 2000)
#merged_Seu <- ScaleData(merged_Seu)
#merged_Seu <- RunPCA(merged_Seu, features = VariableFeatures(object = merged_Seu))
#merged_Seu <- RunUMAP(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- RunTSNE(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- FindNeighbors(merged_Seu, dims = 1:30, verbose = FALSE)
#merged_Seu <- FindClusters(merged_Seu, verbose = FALSE)

#DimPlot(merged_Seu,reduction="umap",group.by="orig.ident",label=TRUE,repel=FALSE)

#FeaturePlot(merged_Seu,reduction="umap",features="mt-Cytb",min.cutoff=0,max.cutoff=4,cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="percent.mt",cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="nFeature_RNA",cols=c("grey","red"))

Save/Load seurat object

#save(merged_Seu, file="merged.Robj")
#load("merged.Robj")

Dataset Integration

#di <- SplitObject(merged_Seu, split.by = "orig.ident")

#for (i in 1:length(di)) {
#  di[[i]] <- NormalizeData(di[[i]], verbose = FALSE)
#  di[[i]] <- FindVariableFeatures(di[[i]], selection.method = "vst", nfeatures = 2000,
#                                      verbose = FALSE)
#}

#dat.anchors <- FindIntegrationAnchors(object.list=di,dims=1:30)
#integrated <- IntegrateData(anchorset=dat.anchors,dim=1:30)

#DefaultAssay(integrated)<-"integrated"
#integrated <- ScaleData(integrated)
#integrated <- RunPCA(integrated,npcs=30)
#integrated <- RunUMAP(integrated,reduction="pca",dims=1:30)
#integrated <- RunTSNE(integrated,reduction="pca",dims=1:30)
#integrated <- FindNeighbors(integrated, dims = 1:30, verbose = FALSE)
#integrated <- FindClusters(integrated, verbose = FALSE)

Load integrated data and create UMAP from original integrated run*

NOTE: loading final object to avoid recalculating ssGSEA data

#load("integrated_orig.Robj")
load("../repo_data/integrated_final.Robj")

Test plots

DefaultAssay(integrated)<-"integrated"
DimPlot(integrated,reduction="umap",split.by="orig.ident",group.by="orig.ident")

DimPlot(integrated,reduction="umap",group.by="orig.ident")

DimPlot(integrated,reduction="umap",group.by="integrated_snn_res.0.8", label=TRUE)

Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
Please use `as_label()` or `as_name()` instead.

Figure 3B - Add cluster/treatment metadata columns and plot labeled UMAP

integrated@meta.data$treat_clust <- paste0(integrated@meta.data$orig.ident,integrated@meta.data$integrated_snn_res.0.8)
integrated@meta.data$clust_treat <- paste0(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
integrated@meta.data$celltype <- integrated@meta.data$integrated_snn_res.0.8
Idents(object = integrated) <- integrated$celltype

integrated <- RenameIdents(object = integrated,  '16' = '16_Tuft')
integrated <- RenameIdents(object = integrated,  '11' = '11_EC')
integrated <- RenameIdents(object = integrated,  '13' = '13_EEC')
integrated <- RenameIdents(object = integrated,  '2' = '2_Stem')
integrated <- RenameIdents(object = integrated,  '5' = '5_Stem')
integrated <- RenameIdents(object = integrated,  '10' = '10_Stem')
integrated <- RenameIdents(object = integrated,  '14' = '14_Paneth')
integrated <- RenameIdents(object = integrated,  '8' = '8_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '9' = '9_Goblet')
integrated <- RenameIdents(object = integrated,  '1' = '1_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '4' = '4_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '15' = '15_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '3' = '3_EC_Progenitor')
integrated <- RenameIdents(object = integrated,  '6' = '6_EC_Progenitor')
integrated <- RenameIdents(object = integrated,  '0' = '0_Early_TA')
integrated <- RenameIdents(object = integrated,  '7' = '7_Early_TA')
integrated <- RenameIdents(object = integrated,  '12' = '12_Unknown')


integrated[["clust_celltype"]] <- Idents(object = integrated)

Fig3b UMAP Plot with color blind safe colors

Idents(object = integrated) <- integrated$clust_celltype

celltype_colors <- c("2_Stem"="#117733",
                     "5_Stem"="#999933",
                     "10_Stem"="#009E73",
                     "0_Early_TA"="#E69F00",
                     "7_Early_TA"="#D55E00",
                     "6_EC_Progenitor"="#0072B2",
                     "3_EC_Progenitor"="#56B4E9",
                     "11_EC"="#88CCEE",
                     "13_EEC"="#6699CC",
                     "1_Secretory_Progenitor"="#661100",
                     "4_Secretory_Progenitor"="#882255",
                     "8_Secretory_Progenitor"="#CC6677",
                     "15_Secretory_Progenitor"="#AA4499",
                     "14_Paneth"="#332288",
                     "16_Tuft"="#000000",
                     "9_Goblet"="#F0E442",
                     "12_Unknown"="#999999")  

dp.cb <- DimPlot(integrated,reduction="umap", cols=celltype_colors, label=TRUE, repel=TRUE, pt.size=2, label.size=6) + NoLegend()
dp.cb

#pdf('Fig3b.cbSafe.pdf',width=14, height=10)
#dp.cb
#dev.off()

Figure 3C - LGR5 Vln plot with color blind safe colors

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
plotOrder <- c("5","2","10","0","1","3","4","6","7","8","9","11","12","13","14","15","16")

vln_colors <- c("2"="#117733",
                "5"="#999933",
                "10"="#009E73",
                "0"="#E69F00",
                "1"="#661100",
                "3"="#56B4E9",
                "4"="#882255",
                "6"="#0072B2",
                "7"="#D55E00",
                "8"="#CC6677",
                "9"="#F0E442",
                "11"="#88CCEE",
                "12"="#999999",
                "13"="#6699CC",
                "14"="#332288",
                "15"="#AA4499",
                "16"="#000000")  

Idents(integrated) <- factor(Idents(integrated), levels= plotOrder)
vl.cb <- VlnPlot(integrated, cols=vln_colors, idents = , features = c("Lgr5"), pt.size = 0.5, slot="data")
vl.cb

#pdf('Fig3c_cbSafe.pdf',width=14, height=8)
#vl.cb
#dev.off()

Figure S3,B - Table of cell counts in each integrated data cluster and sample

p.integrated <- table(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
round(prop.table(p.integrated,2),3)

    
        al     f    rf rf.rapa
  0  0.169 0.191 0.135   0.178
  1  0.096 0.101 0.135   0.119
  2  0.082 0.116 0.104   0.102
  3  0.099 0.091 0.098   0.101
  4  0.125 0.078 0.082   0.061
  5  0.076 0.080 0.084   0.072
  6  0.076 0.074 0.065   0.077
  7  0.056 0.040 0.096   0.055
  8  0.078 0.053 0.034   0.037
  9  0.028 0.046 0.043   0.060
  10 0.034 0.032 0.044   0.055
  11 0.026 0.056 0.033   0.029
  12 0.039 0.028 0.024   0.032
  13 0.005 0.006 0.009   0.009
  14 0.005 0.005 0.006   0.007
  15 0.002 0.002 0.005   0.005
  16 0.004 0.002 0.004   0.001

Figure 5D - heatmap

DefaultAssay(integrated)<- "integrated"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
dd <- subset(integrated, idents = c("2", "5", "10"), downsample=100)

topvst <- head(VariableFeatures(dd), 500)
mat <- as.matrix(dd@assays$integrated@scale.data) #as.matrix(subset_dd@assays$integrated@scale.data)
mat <- mat[topvst,]

genes <- c(GOI.lists$Biton_S1_ISC.I, GOI.lists$Biton_S1_ISC.II, GOI.lists$Biton_S1_ISC.III)
labels <- c(rep('Biton ISC I', length(GOI.lists$Biton_S1_ISC.I)), 
            rep('Biton ISC II', length(GOI.lists$Biton_S1_ISC.II)), 
            rep('Biton ISC III', length(GOI.lists$Biton_S1_ISC.III)))

idxs <- which(genes %in% rownames(mat))
genes <- genes[idxs]
labels <- labels[idxs]
mat <- mat[genes,]

ht <- Heatmap(mat, column_names_side = 'top', row_names_gp = gpar(fontsize = 9), column_names_gp = gpar(fontsize = 12),
              column_title = '', name = 'scaled data', column_labels = rep('', ncol(mat)),
              column_split = factor(as.character(dd$integrated_snn_res.0.8), levels = c('5', '2', '10')), 
              row_split = labels, row_order = genes, #column_order = sort(colnames(mat)),
              cluster_column_slices = FALSE,
              top_annotation = HeatmapAnnotation(cluster = as.character(dd$integrated_snn_res.0.8)))

draw(ht)

#pdf('Fig3D.pdf',width=12, height=10)
#draw(ht)
#dev.off()

GSEA dot plots

data <- as.data.frame(read_excel('cluster_2_5_10stem_gsea.xlsx', sheet = 'Sheet2'))
data

data$FDR <- data$`FDR q-val` + 0.001
data$Comparison <- gsub('\\.noNA\\.NES','',data$Comparison)

comparisons <- c('f5_v_al5', 'rf5_v_al5', 'rf.rapa5_v_al5')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl5 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) + 
  geom_point() + 
  scale_color_gradient2(midpoint=0, low="blue", mid="white",
                        high="red", space ="Lab", limits=c(-3,3))+
  scale_size_continuous(range=c(1,6))+
  ggtitle('Cluster 5 GSEA') + theme_classic() + 
  theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl5

#pdf('Fig5F.pdf',width=4, height=4)
#cl5
#dev.off()

comparisons <- c('f2_v_al2', 'rf2_v_al2', 'rf.rapa2_v_al2')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl2 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) + 
  geom_point() + 
  scale_color_gradient2(midpoint=0, low="blue", mid="white",
                        high="red", space ="Lab", limits=c(-3,3))+
  scale_size_continuous(range=c(1,6))+
  ggtitle('Cluster 2 GSEA') + theme_classic() + 
  theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl2

#pdf('FigS3Da.pdf',width=4, height=4)
#cl2
#dev.off()

comparisons <- c('f10_v_al10', 'rf10_v_al10', 'rf.rapa10_v_al10')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl10 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) + 
  geom_point() + 
  scale_color_gradient2(midpoint=0, low="blue", mid="white",
                        high="red", space ="Lab", limits=c(-3,3))+
  scale_size_continuous(range=c(1,6))+
  ggtitle('Cluster 10 GSEA') + theme_classic() + 
  theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl10

#pdf('FigS3Db.pdf',width=4, height=4)
#cl10
#dev.off()

#there is a glitch in this plot, cl2 loses x axis, legend order is different

figS3D <- ggarrange(cl2,cl10, ncol=2, nrow=1)
figS3D

#pdf('FigS3D_incorrect.pdf',width=8, height=4)
#figS3D
#dev.off()

Extended data Figure S4A - Feature Plots - color blind safe

DefaultAssay(integrated)<-"RNA"
colorScheme <- c("#C1BEBF","#fe6100")

fp.Muc2 <- FeaturePlot(integrated,reduction="umap",features="Muc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Muc2

fp.Tff3 <- FeaturePlot(integrated,reduction="umap",features="Tff3",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Tff3

fp.Lyz1 <- FeaturePlot(integrated,reduction="umap",features="Lyz1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lyz1

fp.Defa30 <- FeaturePlot(integrated,reduction="umap",features="Defa30",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Defa30

fp.Chga <- FeaturePlot(integrated,reduction="umap",features="Chga",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Chga

fp.Reg3a <- FeaturePlot(integrated,reduction="umap",features="Reg3a",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Reg3a

fp.Alpi <- FeaturePlot(integrated,reduction="umap",features="Alpi",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Alpi

fp.Atoh1 <- FeaturePlot(integrated,reduction="umap",features="Atoh1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Atoh1

fp.Lgr5 <- FeaturePlot(integrated,reduction="umap",features="Lgr5",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lgr5

fp.Smoc2 <- FeaturePlot(integrated,reduction="umap",features="Smoc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Smoc2

figS4A <- ggarrange(fp.Muc2,fp.Tff3,fp.Lyz1,fp.Defa30,fp.Chga,fp.Reg3a,fp.Alpi,fp.Atoh1,fp.Lgr5,fp.Smoc2, legend = "right", ncol = 5, nrow = 2)

figS4A

#pdf('FigS4A_cbSafe.pdf',width=20, height=10)
#figS4A
#dev.off()

Figure 5E - cell cycle data - color blind safe

DefaultAssay(integrated)<-"RNA"
gene.names <- rownames(integrated@assays$RNA@data)
Idents(object = integrated) <- integrated$seurat_clusters
stem.subset <- subset(integrated, idents = c("2","5","10"))
levels(stem.subset) <- c("5","2","10")

vln_colors <- c("2"="#117733",
                "5"="#999933",
                "10"="#009E73")
  
s.InData <- intersect(gene.names,GOI.lists$mm.s)
g2m.InData <- intersect(gene.names,GOI.lists$mm.g2m)
stem.subset[["percent.mm.s"]] <- PercentageFeatureSet(stem.subset, features = s.InData)
stem.subset[["percent.mm.g2m"]] <- PercentageFeatureSet(stem.subset, features = g2m.InData)
mm.s <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.s",pt.size = 0.3,slot = "data")
mm.s

mm.g2m <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.g2m",pt.size = 0.3,slot = "data")
mm.g2m

figS4b <- ggarrange(mm.s,mm.g2m, legend = FALSE, ncol=2, nrow=1)
figS4b

pdf('figS4b_cbSafe.pdf',width=8, height=4)

figS4b
dev.off()

png 
  2 

integrated[["percent.mm.s"]] <- PercentageFeatureSet(integrated, features = s.InData)
integrated[["percent.mm.g2m"]] <- PercentageFeatureSet(integrated, features = g2m.InData)

mm.g2m.Flag <- if_else(integrated@meta.data$percent.mm.g2m >= 0.3, "Yes", "No")
integrated@meta.data$mm.g2m.Flag <- mm.g2m.Flag

p <- table(integrated$clust_treat,integrated$mm.g2m.Flag)
p.g2m.summary <- round(prop.table(p,2),3)

mm.s.Flag <- if_else(integrated@meta.data$percent.mm.s >= 0.2, "Yes", "No")
integrated@meta.data$mm.s.Flag <- mm.s.Flag

p <- table(integrated$clust_treat,integrated$mm.s.Flag)
p.s.summary <- round(prop.table(p,1),3)

escape ssGSEA - run ssGSEA to quantify expression of the BitonI gene set clusters

This code is no longer working due to R and package updates but resulting data is stored in the seurat object escape 1.0.0 is probably required but is no longer available. the version in this R image is 1.0.1 escape ssGSEA - run ssGSEA to quantify expression of the BitonI gene set clusters

egs <- GeneSet(GOI.lists$Biton_S1_ISC.I, setName="BitonI")
ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)

[1] "Using sets of 1000 cells. Running 19 times."

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘gsva’ for signature ‘"matrix", "character"’

escape ssGSEA - run ssGSEA to quantify expression of the BitonII gene set clusters

# egs <- GeneSet(GOI.lists$Biton_S1_ISC.II, setName="BitonII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonII_ssGSEA <- ES$BitonII

escape ssGSEA - run ssGSEA to quantify expression of the BitonIII gene set clusters

# egs <- GeneSet(GOI.lists$Biton_S1_ISC.III, setName="BitonIII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonIII_ssGSEA <- ES$BitonIII

Extended data 5 plots - colorblind safe

Figure 5G - Biton 1 in cluster5 Figure 5G_II - Biton II in cluster2 Figure 5G_III - Biton III in cluster10 Figure 5H - Pdgfa in cluster5 Figure S3E - Gkn3 in cluster5*

Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust

vln_colors <- c("al5"="#e69f00",
                "f5"="#56b4f9",
                "rf5"="#117733",
                "rf.rapa5"="#d55e00")

plot.order <- c("al5","f5","rf5","rf.rapa5")
Idents(object = clust5.subset) <- factor(Idents(object = clust5.subset), levels = plot.order)
vln.BitonI.ssGSEA <- VlnPlot(clust5.subset, cols = vln_colors, features="BitonI_ssGSEA",slot = "data",pt.size = 0.4) + 
  stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-I_ssGSEA")

vln.BitonI.ssGSEA

#pdf('Fig5G.pdf',width=4, height=4)
#vln.BitonI.ssGSEA
#dev.off()

vln.Pdgfa <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Pdgfa"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Pdgfa

pdf('FigS4g_cbSafe.pdf',width=4, height=4)

vln.Pdgfa
dev.off()

png 
  2 

vln.Gkn3 <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Gkn3"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Gkn3

pdf('FigS4h_cbSafe.pdf',width=4, height=4)

vln.Gkn3
dev.off()

png 
  2 

clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust

vln_colors <- c("al2"="#e69f00",
                "f2"="#56b4f9",
                "rf2"="#117733",
                "rf.rapa2"="#d55e00")

plot.order <- c("al2","f2","rf2","rf.rapa2")
Idents(object = clust2.subset) <- factor(Idents(object = clust2.subset), levels = plot.order)
vln.BitonII.ssGSEA <- VlnPlot(clust2.subset, cols = vln_colors, features="BitonII_ssGSEA",slot = "data",pt.size = 0.4) + 
  stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-II_ssGSEA")
vln.BitonII.ssGSEA

#pdf('Fig5G_II.pdf',width=4, height=4)
#vln.BitonII.ssGSEA
#dev.off()

clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust

vln_colors <- c("al10"="#e69f00",
                "f10"="#56b4f9",
                "rf10"="#117733",
                "rf.rapa10"="#d55e00")

plot.order <- c("al10","f10","rf10","rf.rapa10")
Idents(object = clust10.subset) <- factor(Idents(object = clust10.subset), levels = plot.order)
vln.BitonIII.ssGSEA <- VlnPlot(clust10.subset, cols = vln_colors, features="BitonIII_ssGSEA",slot = "data",pt.size = 0.4) + 
  stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-III_ssGSEA")
vln.BitonIII.ssGSEA

#pdf('Fig5G_III.pdf',width=4, height=4)
#vln.BitonIII.ssGSEA
#dev.off()

Biton plots assembled

bitonPlot <- ggarrange(vln.BitonI.ssGSEA,vln.BitonII.ssGSEA,vln.BitonIII.ssGSEA,nrow=1,ncol=3)
bitonPlot 

#pdf('FigS4f_cbSafe.pdf',width=10, height=5)
#bitonPlot
#dev.off()

Fig 5H style plots for Oat

Figure 5H style plots for Oat, Oct, Ass Asl in cluster5

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust
levels(clust5.subset) <- c("al5","f5","rf5","rf.rapa5")
clust5.subset$treat_clust <- factor(x = clust5.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5"))
VlnPlot(clust5.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

The default behaviour of split.by has changed.
Separate violin plots are now plotted side-by-side.
To restore the old behaviour of a single split violin,
set split.plot = TRUE.
      
This message will be shown once per session.

vln.Oat.5 <- VlnPlot(clust5.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5

#pdf('Fig5H_Oat.5.pdf',width=4, height=4)
#vln.Oat.5
#dev.off()

clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust
levels(clust2.subset) <- c("al2","f2","rf2","rf.rapa2")
clust2.subset$treat_clust <- factor(x = clust2.subset$treat_clust, levels = c("al2","f2","rf2","rf.rapa2"))
VlnPlot(clust2.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

vln.Oat.2 <- VlnPlot(clust2.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.2

#pdf('Fig5H_Oat.2.pdf',width=4, height=4)
#vln.Oat.2
#dev.off()

clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust
levels(clust10.subset) <- c("al10","f10","rf10","rf.rapa10")
clust10.subset$treat_clust <- factor(x = clust10.subset$treat_clust, levels = c("al10","f10","rf10","rf.rapa10"))
VlnPlot(clust10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

#vln.Oat.10 <- VlnPlot(clust10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
#vln.Oat.10

#pdf('Fig5H_Oat.10.pdf',width=4, height=4)
#vln.Oat.10
#dev.off()

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.2.10.subset <- subset(integrated, idents = c("5","2","10"))
Idents(object = clust5.2.10.subset) <- clust5.2.10.subset$treat_clust
levels(clust5.2.10.subset) <- c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10")
clust5.2.10.subset$treat_clust <- factor(x = clust5.2.10.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10"))
VlnPlot(clust5.2.10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

vln.Oat.5.2.10 <- VlnPlot(clust5.2.10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5.2.10

#pdf('Fig5H_Oat.5.2.10.pdf',width=12, height=4)
#vln.Oat.5.2.10
#dev.off()

Figure 5H style plots but without cluster for Oat*

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$orig.ident
levels(integrated) <- c("al","f","rf","rf.rapa")
integrated$orig.ident <- factor(x = integrated$orig.ident, levels = c("al","f","rf","rf.rapa"))
VlnPlot(integrated, features = c("Oat"), split.by = "orig.ident", group.by = "orig.ident", slot= "data",pt.size = 0.5)

vln.Oat.noclust <- VlnPlot(integrated, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.noclust

#pdf('Fig5H_Oat.noclust.pdf',width=4, height=4)
#vln.Oat.noclust
#dev.off()

Figure S3E - Differential Expression Testing to prepare .rnk files*

# DefaultAssay(integrated)<-"RNA"
# Idents(object = integrated) <- integrated$treat_clust
# levels(x = integrated)
# ##
# #adjust the id1 and id2 variables to set up different tests
# #could use a loop for this
# ##
# id1 <- "al5"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
# 
# id1 <- "al10"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
# 
# id1 <- "al10"
# id2 <- "al5"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)

Save integrated object with ssGSEA data

#save(integrated, file="integrated_final.Robj")

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8      
 [8] LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ComplexHeatmap_2.6.2        scater_1.18.6               SingleCellExperiment_1.12.0 GSEABase_1.52.1             graph_1.68.0                annotate_1.68.0            
 [7] XML_3.99-0.5                AnnotationDbi_1.52.0        dittoSeq_1.2.6              escape_1.0.1                DESeq2_1.30.1               SummarizedExperiment_1.20.0
[13] Biobase_2.50.0              MatrixGenerics_1.2.1        matrixStats_0.58.0          GenomicRanges_1.42.0        GenomeInfoDb_1.26.7         IRanges_2.24.1             
[19] S4Vectors_0.28.1            BiocGenerics_0.36.1         kableExtra_1.3.1            knitr_1.31                  sctransform_0.3.2           ggpubr_0.4.0               
[25] ggrepel_0.9.1               ggplot2_3.3.3               openxlsx_4.2.3              readxl_1.3.1                Matrix_1.2-18               dplyr_1.0.4                
[31] umap_0.2.7.0                cowplot_1.1.1               Seurat_3.2.3               

loaded via a namespace (and not attached):
  [1] reticulate_1.18           tidyselect_1.1.0          RSQLite_2.2.3             htmlwidgets_1.5.3         BiocParallel_1.24.1       Rtsne_0.15                munsell_0.5.0            
  [8] codetools_0.2-16          ica_1.0-2                 future_1.21.0             miniUI_0.1.1.1            withr_3.0.0               colorspace_2.0-0          rstudioapi_0.13          
 [15] ROCR_1.0-11               ggsignif_0.6.0            tensor_1.5                listenv_0.8.0             labeling_0.4.2            GenomeInfoDbData_1.2.4    polyclip_1.10-0          
 [22] farver_2.0.3              pheatmap_1.0.12           bit64_4.0.5               parallelly_1.23.0         vctrs_0.6.5               generics_0.1.0            xfun_0.21                
 [29] R6_2.5.1                  clue_0.3-58               ggbeeswarm_0.6.0          rsvd_1.0.3                msigdbr_7.2.1             locfit_1.5-9.4            bitops_1.0-6             
 [36] spatstat.utils_2.0-0      cachem_1.0.3              DelayedArray_0.16.3       assertthat_0.2.1          promises_1.1.1            scales_1.1.1              beeswarm_0.2.3           
 [43] gtable_0.3.0              beachmat_2.6.4            Cairo_1.5-12.2            globals_0.14.0            goftest_1.2-2             rlang_1.1.4               genefilter_1.72.1        
 [50] GlobalOptions_0.1.2       splines_4.0.3             rstatix_0.6.0             lazyeval_0.2.2            broom_0.7.4               yaml_2.2.1                reshape2_1.4.4           
 [57] abind_1.4-5               backports_1.2.1           httpuv_1.5.5              tools_4.0.3               ellipsis_0.3.2            RColorBrewer_1.1-2        ggridges_0.5.3           
 [64] Rcpp_1.0.6                plyr_1.8.6                sparseMatrixStats_1.2.1   zlibbioc_1.36.0           purrr_1.0.2               RCurl_1.98-1.2            rpart_4.1-15             
 [71] openssl_2.2.0             deldir_0.2-9              GetoptLong_1.0.5          viridis_0.5.1             pbapply_1.4-3             zoo_1.8-8                 haven_2.3.1              
 [78] cluster_2.1.0             magrittr_2.0.3            data.table_1.13.6         RSpectra_0.16-0           scattermore_0.7           circlize_0.4.12           lmtest_0.9-38            
 [85] RANN_2.6.1                fitdistrplus_1.1-3        GSVA_1.38.2               hms_1.0.0                 patchwork_1.1.1           mime_0.9                  evaluate_0.24.0          
 [92] xtable_1.8-4              rio_0.5.16                shape_1.4.5               gridExtra_2.3             compiler_4.0.3            tibble_3.0.6              KernSmooth_2.23-17       
 [99] crayon_1.4.1              htmltools_0.5.8.1         mgcv_1.8-33               later_1.1.0.1             tidyr_1.1.2               geneplotter_1.68.0        DBI_1.1.1                
[106] MASS_7.3-53               car_3.0-10                cli_3.6.3                 igraph_1.2.6              forcats_0.5.1             pkgconfig_2.0.3           foreign_0.8-80           
[113] scuttle_1.0.4             plotly_4.9.3              xml2_1.3.2                vipor_0.4.5               webshot_0.5.2             XVector_0.30.0            rvest_0.3.6              
[120] stringr_1.4.0             digest_0.6.36             RcppAnnoy_0.0.18          spatstat.data_2.0-0       rmarkdown_2.6             cellranger_1.1.0          leiden_0.3.7             
[127] edgeR_3.32.1              uwot_0.1.10               DelayedMatrixStats_1.12.3 curl_5.2.1                shiny_1.6.0               rjson_0.2.20              lifecycle_1.0.4          
[134] nlme_3.1-149              jsonlite_1.8.8            BiocNeighbors_1.8.2       carData_3.0-4             limma_3.46.0              viridisLite_0.3.0         askpass_1.1              
[141] pillar_1.4.7              lattice_0.20-41           fastmap_1.2.0             httr_1.4.2                survival_3.2-7            glue_1.4.2                zip_2.1.1                
[148] spatstat_1.64-1           png_0.1-7                 bit_4.0.4                 stringi_1.5.3             blob_1.2.1                BiocSingular_1.6.0        memoise_2.0.1            
[155] irlba_2.3.3               future.apply_1.7.0       

writeLines(capture.output(sessionInfo()), "2024_sessionInfo.txt")